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Abstract 

We present a theoretical framework using quorum-percolation for describing the initiation 
of activity in a neural culture. The cultures are modeled as random graphs, whose nodes are 
excitatory neurons with k- m inputs and fc out outputs, and whose input degrees fcj n = k obey 
given distribution functions pk- We examine the firing activity of the population of neurons 
according to their input degree (fc) classes and calculate for each class its firing probability 
<3>fc(i) as a function of t. The probability of a node to fire is found to be determined by its 
in-degree fc, and the first-to-fire neurons are those that have a high k. A small minority of 
high-fc classes may be called "Leaders," as they form an inter-connected subnetwork that 
consistently fires much before the rest of the culture. Once initiated, the activity spreads 
from the Leaders to the less connected majority of the culture. We then use the distribution 
of in-degree of the Leaders to study the growth rate of the number of neurons active in 
a burst, which was experimentally measured to be initially exponential. We find that this 
kind of growth rate is best described by a population that has an in-degree distribution 
that is a Gaussian centered around k — 75 with width a = 31 for the majority of the 
neurons, but also has a power law tail with exponent —2 for ten percent of the population. 
Neurons in the tail may have as many as k = 4, 700 inputs. We explore and discuss 
the correspondence between the degree distribution and a dynamic neuronal threshold, 
showing that from the functional point of view, structure and elementary dynamics are 
interchangeable. We discuss possible geometric origins of this distribution, and comment 
on the importance of size, or of having a large number of neurons, in the culture. 

Keywords: Neuronal cultures, Graph theory, Activation dynamics, Percolation, Statistical 
mechanics of networks, Leaders of activity, Quorum. 
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I. INTRODUCTION 



Development of connectivity in a neuronal network is strongly dependent upon the environment 
in which the network grows: cultures grown in a dish will develop very differently from networks 
formed in the brain. In a dish, the only signals that neurons are exposed to are chemicals secreted 
by neighboring neurons, which then must diffuse to other neurons via the large volume of fluid that 
surrounds the culture. The result is a connectivity dominated by proximity in a pl anar geometry , 
whos e input degree follows a statistical distribution function that is Gaussian-like (|Soriano et al. . 



2008[). This is contrasted by the intricate guidance of axons during the creation of connectivity 



in the brain, which is dictated by a detailed and very complex "blueprint" for connectivity. As a 
result, the firing pattern of a culture is an all-or-none event, a population spike in which practically 
all the neurons participate and are simultaneously active for a brief period of time, spiking about 
3-4 times on average. 



We have previously shown (ICohen et a/1 1201 0|) that graph theory and statistical mechanics are 



useful in unraveling properties of the network in a rat hippocampal culture, mostly because of the 
statistical nature of the connectivity. With these tools we have been able to understand such phe- 
nomena as the degree distribution of in put connections, the ratio of inhibitory to e xcitatory neurons 



2006; 



Soriano et all 



2008). We found that a 



and the input cluster size distribution (IBreskin et all 
Gaussian degree distribution gives a good quantitative description for statistical properties of the 
network such as the appearance of the giant m-connected component and its size as a function of 
connectivity. The inhibitory component was found to be about 30% in hippocampal cultures, and 
about 20% in cortex. We furthermore obs erved that the appea rance of a fully connected network 



coincides precisely with the time of birth (|Soriano et al. 



2008) 



In this paper we apply our graph theoretic approach to the intriguing process of the initiation of 
a spontaneous population spike. On the one hand, a perturbation needs to be created that pushes a 
number of neurons to begin firing. On the other hand, the initial firing pattern must propagate to 
the rest of the neurons in the culture. Understanding this recruitment process will give insight on 
the structure of the network, on the interrelation of activation in neurons, and on the dynamics of 
neuronal firing in such a complex culture. A simple scenario for initiation that one might conceive 
of is wave-front propagation, in which a localized and limited area of initiation is ignited first, and 
from there sets up a spherical traveling front of excitation. However, as we shall see, the initiation 
is a more intricate process. 
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The experimental situation regarding initiation of activity is complex. In quasi one-dimensional 
networks we have been able to show that activity originates in a single "Burst Initiation 
Zone" (BIZ), in which a slow recruitment process occurs oyer several hundreds of millis econds 



(IFeinerman et all 



2005 



Golomb and Ermentrout, 



1999, 



2002; 



Osan and Ermentrout! 



2002). From 



this BIZ the activity propagates to the res t of the linear cu l ture a long an orderly and causal path 



dictated by the one dimensional structure(Feinerman et al. 



such causal propagation is not observed (IMaeda et al 



1995 



2007). In two dimensional networks 



Streit et al. 



20011). and the precise 



mode of propagation has not been identified. Recently, Eytan and Marom (lEytan and Marom , 



2006|) found that a small subgroup of neurons were the "first-to-fire," and that also in this case the 



initiation is long (on the order of hundreds of milliseconds). They also observed that the growth 
rate is exponential in the initial stage and then changes to faster than exponen tial. These neurons 



were later shown to characterize and "lead" the burst (Eckmann et al 



2008|) . and to recruit the 



2008). 



neurons in their proximity in the "pre-burst" period (E ckmann et all 

In this paper we address and connect three experimental observations that are at first sight 
unrelated. The first and fundamental observati on is the fact that bursts are initiated by Leaders, or 
first-to-fire neurons ( Eytan and Marom . I2OO6 ). We use the quorum-percolation model to answer 
the question - what makes these Leaders different from the rest of the network - by showing that 
one of their important characteristics is a high in-degree, i.e., a large number of input connections. 

We then turn to the second experimental observation, that the activity in a burst starts with 
an exponential growth. We show that this can happen only if the distribution of in-degree is a 
power law with exponent of —2. However, this needs to be related with the third experimental 
observation, which is that the distribution of in-degrees is Gaussian rather than power law. We 
reconcile both observations by stitching together the two solutions into an in-degree distribution 
that has the large majority of neurons in a Gaussian centered around an average in-degree of 
about 75, while ten percent of the population lie on a power law tail that can reach a few thousand 
connections. We show that this reconstructs the full experimentally observed burst structure, which 
is an exponential initiation during the pre-burst followed by a super-exponential during the burst. 

We present several ideas on the origin of these distributions in the spatial extent and geometry 
of the neurons, and show that the distribution of in-degrees is proportional to the distribution of 
spatial sizes of the dendritic trees. We thus conjecture that the distribution of dendritic trees is 
mostly Gaussian, but that a few neurons must have dendrites that go off very far, with power law 
distribution of this tail. We end by making some additional conjectures about one-dimensional 
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cultures and on the importance of the size of the culture. 



II. METHODS 



A. Quorum percolation model for dynamics of random graph network 



We describe the neural culture using a simplified model of a network whose nodes are neurons 
with links that are the axon/dendri te connections. T his picture is further simplified if we consider a 
randomly connected sparse graph (IBollobasl . 1 1998b . with a uniform strength on all the connections. 
The structure and topology of the graph are determined by specifying a probability distribution p k 
for a node to have k inputs. Percolation on a graph is the process by which a property spreads 
through a sizeable fraction of the graph. In our case, this property is the firing of neurons. The 
additional characteristic of Quorum Percolation is, as its name implies, that a burst of firing activity 
will propagate throughout the neuronal culture only if a quorum of more than m firing neighbors 
has ignited on the corresponding graph. While this description makes a number of assumptions 
that are not exact in their comparison to the experimen t, it does, as we have been able to show 



previ ously, capture the essential behavior of the network (ICohen et all 



20101 : 



Tlustv and Eckmann , 



2009J). The use of such a simple model for the neuronal network is justified at the end of the 
Methods section. 

In particular, in this paper we obtain a theoretical explanation for the experimental observation 
of initiation of activity by a small number of neurons from the network and the subsequent gradual 
recruitment of the rest of the network. Within the framework of a random graph description, we 
have previously shown that the dynamics of firing in the network is described by a fixed point 
equation for the probability of firing in the network, which also corresponds to the experimentally 



observed fractio n of neurons that fire . Experimenta 



external voltage (IBreskin et al. . 



2006; 



Soriano et al. . 



ly, th is fraction can be set by applying an 



2008). In the case of spontaneous activity, 



this fraction is determined by the interplay of noise and the intrinsic sensitivity of the neurons 



(lAlvarez-Lacalle and Moses . 



2009) 



Specifically, connections are described by the adjacency matrix A, produced according to the 
probability distribution p k , with = 1 if there is a directed link from j to i, and A^ = 
otherwise and An = 0. 

Our study starts by assuming initial conditions where a fraction / of neurons are switched "on" 
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externally at time t — 0. The neurons fire, and once they do they stay on forever - a neuron will 
be on at time t + 1 if at time t it was on. A neuron is either turned on at t = (with probability 
/) or, if it is off at time t, then it will it will be turned on at time t + 1 if at least m of its upstream 
(incoming) nodes were on at time t: 

Si (t + 1) = Si (t) + (1 - Si (t)) 6 A H s j(t) - m j , (!) 

where Sj(t) describes the state of the neuron at time t (1 for "on" and for "off") and 6 is the 
step function (1 for x > 0, otherwise). Note that the second term, which accounts for the 
neuron's probability to fire at time step t, creates a coupling of Si(t) to all its inputs Sj(t). Lacking 
a turning off process, the number of active nodes is monotonically increasing and converges into 
a steady-state (a fixed point) within a finite time tf, which is smaller than the number of neurons, 
t f <N. 

To derive the "mean-field" dynamics for the average fraction of firing neurons, <3>(i) = 
(si(t)) = A r_1 J2i one cannot simply average equation (OQ) directly. This is due to the corre- 
lation between Sj(t) and AijSj(t) — rnj . The correlation exists because if a given neuron 
fires, Si(t) = 1, and it was not externally excited, then at least m of its inputs are firing and the 
step function over its inputs must also be 1. In fact, at the fixed point the correlation is strong, 
(1 — Siit)) 6 (^2j AijSj{t) — rnj = 0, since in the steady state a neuron can remain off if and 
only if less than m of its inputs fire. To avoid the correlations, one utilizes the monotonicity to 
realize that a neuron is on only if it was turned on externally at time t = or, if it was off at 
t = then at some time t at least m of its inputs fired. We can then replace Si(t) by Sj(0) and 
rewrite equation CO) as Sj(t + 1) = Sj(0) + (1 — Sj(0)) 6 (j^j AijSj(t) — rnj . In the tree approx- 
imation, which disregards loop feedbacks, the initial firing state of a neuron, Sj(0), cannot affect 
its inputs, Sj(t). Inversely, it is obvious that Sj(0) which is determined externally is independent 
of Sj(t). Therefore, Sj(0) and (j^j AijSj(t) — raj can be averaged independently. The result is 
the mean- field iteration map 

$(* + l) = /+(l-/)¥(m,$(t)), (2) 

where the combinatorial expression for \1/ is the probability that at least m inputs are firing and / 
is the initial firing fraction, / = $(0) = (sj(0)). The steady state of the network is defined by the 
fixed point $oo, which is found by inserting $(£ + 1) = $(£) = into equation ©, to obtain 

$oo = /+(l-/)*(m,$ 00 ). 
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B. Collectivity and the critical point 



The role of Leader neurons in the initiation and the development of bursts can be clarified by 
dividing the neurons into classes of in-degree fc ("fc-class") and looking at the dynamics of each 
class separately. The total firing probability $ = ^ k Pk^k is thus composed of the sum over the 
individual probabilities $^ for each fc-class to fire. The mean field equation for a given fc-class is 



where ^ k is the probability that a neuron with fc inputs has at least m that are firing. Although all 
fc-classes are coupled through the common $, the formulation of equation [3] allows the tracing of 
the fraction of each class and its dynamics during its evolution. 

It follows from © that at the fixed point $fc >00 = / + (1 — /)W fc (m, $oo)- Given the time 
dependence of one can extract the fraction of firing neurons in each fc-class, by plugging 
into ©. 

The combinatorial expressions for ^ and ^ are: 



There is a particular critical initial firing /* where the solution jumps from $ ~ / (i.e., practi- 
cally all activation is externally driven and there is almost no collectivity, ^ <C 1) to $ ~ 1 where 
most firing is due to inputs (and ^ ~ 1). It is both convenient and instructive to treat and simulate 
the network near this transition point, since the dynamics there is slow. This allows the different 
steps in the recruitment process to be easily identified and distinguished. 

C. Simulation of the quorum percolation model 

The model was numerically investigated by performing by a simulation, employing N = 
500, 000 neurons. This number was chosen to match as close as possible the number of neu- 
rons typically in an experiment, which is on the order of one million. An initial fraction / of the 
neurons were randomly selected and set to "on" (i.e., fired). At every time step a neuron would 
fire if it fired the previous step, or if more than a threshold number of its input neighbors fired. 
The threshold m and the initial firing component / could be varied, and the activity history of all 
neurons was stored for subsequent analysis. We used a number of different degree distributions, 



$ fc (* + l) = /+(l-/)¥ fc (m, $(*)), 



(3) 




(4) 
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including a Gaussian, exponential and power law for the network. We also used a tailored Gaus- 
sian distribution in which 10% of the high k neurons, which have k higher than a given k tai \, obey 
a power law distribution function. 



D. Validity of the model 

1. The use of a random graph for neuronal cultures 

The spatial extent and arrangement of connections can be of importance to the dynamics of the 
network. In contrast to spatially embedded (metric) graphs, random graphs allow any two nodes 
to connect, z'.e.they are the analog of infinite dimensional networks. The experiment is obviously 
metric but our model employs a random g raph. This seeming contradiction was resolved in a 



previous study (ITlusty and Eckmannl . 120091) . where we showed that if the average connectivity is 
high enough then the graph is effectively random (i.e. of very high dimension). Why does the 
random graph picture describe so successfully the measurements of a 2D neural culture while it 



comp letely neglects the notions of space and vicinity? As we explained in (|Tlusty and Eckmann , 



20091) . there is a basic difference in the manner in which random and metric graphs are ignited. 
In metric graphs it suffices to initially turn on localized excitation nuclei, which are then able to 
spread an excitation front throughout the spatially extended network. In random graphs, there are 
no such nuclei an d one has to excite a finite fraction of the neurons to keep the ignition going. 



Still, we showed (ITlusty and Eckmann , 



20091) that the experimental network, which is obviously 
an example of a metric graph, is effectively random, since its finite size and the demand for a large 
quorum of firing inputs makes the occurrence of excitation nuclei very improbable. As explained 
in that paper, this occurs when it becomes impossible to identify causal paths in space along which 
activity propagates, with one neuron activating its neighbor and so on. In other words, the activity 
burst does not initiate at one specific nucleating site, and has instead multiple locations at which 
activity appears. This is exactly the characteristic of a random graph with no spatial correlations. 
Thus a highly connected graph in 2D such as ours has characteristics that are similar to a high 
dimensional graph with near-neighbor connectivity. 
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2. Approximating neuronal cultures as tree-like graphs 



The basic reason why a tree-like graph will describe the experimental network lies in the ob- 
servation that the percentage of connections eman ating from a neuron t hat actually participate in 
a loop is small. Indeed, we have demonstrated in (ISoriano et all 120081) that the average number 
of connections per neuron is large, about 100. On the other hand, because the network is built 
from dissociated neurons, the connectivity is determined by a spatial search process during their 
growth, which is for all practical purposes a random one. We therefore have a random network 
(embedded in a metric space) with about 100 connections per node. 

Such networks do indeed have some loops, and thus we need to study the effect of such loops 
on the general argument of equation [2] For this, it suffices to study the case of 2-loops (which in 
fact cause the strongest correlations). Assume neuron A and neuron B are linked in a 2-loop. 

If neuron A fires at time t then it does not change any more, and therefore the state of B at time 
t + 1 does not matter for the state of A at any later time. If A is not firing at time t then it decreases 
the probability of B to fire at time t + 1, and this in turn reduces its own probability to fire at time 
t + 2. Clearly, this effectively decreases the probability of A to fire. We shall now show that this 
effect is 1/ [k 2 N) where k is the mean degree (100) and N is the total number of neurons that B 
can connect to. 

To see this, we note that if A is off then the number of available inputs that can fire into B 
reduces by one, from k to k — 1. We show below that this corresponds changing the ignition level 
$ from m/k to m/(k — 1). The overall effect of a 2-loop on the probability of B to fire therefore 
scales like m/(k — 1) — m/k ~ m/k 2 . The back-propagated effect on A will be of order m 2 /k A . 

To estimate the total number of 2-loops that include neuron A we first look at all trajectories of 
length 2 that emanate from A. There are k 2 nt such trajectories. Of these a fraction of k in /N will 
return to A. The total number of 2-loops that start at A is thus k 2 ut k^/N . The fraction of inputs of 
A that participate in a 2-loop is therefore k 2 ut /N 

Assuming that for a typical neuron k in is equal to A; out , the total effect of 2-loops on $ calculated 
at neuron A is therefore m 2 k 2 /N/ k 4 = m 2 / (k 2 N). 

Below we show the applicability of the tree like random graph model by comparing its results 
directly to the simulation that uses iV = 500, 000 neurons. The correspondence between model 
and simulations (Figures 0] and [5]) is satisfactory. 

In the experimental case, spatial proximity may lead to more connections than in a random 
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graph. The effect of space is to change the relevant number of neurons N from the total number 
to those that are actually accessible in 2D. That number N is on the order of N spiice = 3, 000 as 
compared to iV tota i = 500, 000. However, l/k 2 N is still a small number. 



2010bl) . 



In a separate work, a simulation that takes space into account was performed (|ZbindenL 
and the number of loops could be evaluated directly. Indeed we found that the number of loops is 
enhanced over the random graph estimation by a factor of A^ pacc /A/t tai» but still remains small. 



3. Applicability of the averaged equation (mean-field) approach 



In a physical model one must be sure that the ensemble of random examples chosen to average 
over a given quantity does indeed represent well the statistics of the system that is being treated. In 
our model, the connectivity of the graph is fixed ("quenched" disorder), and the ensemble is that of 
the random graphs that can be generated with the particular choice of input connection distribution 
function. In reality, the experiment and the simulation measure the bursts inside one particular 
realization. However, the mean-field equation averages over a whole ensemble of such random 
graphs. The question is whether the averages obtained using one real graph are representative of 
the whole ensemble. This is a behavior which is termed "self-averaging", and means that, in the 
limit of large graphs, one single configuration represents the average behavior of the ensemble. 

The similar ensemble of the classic Hopfield (spin-gla ss) model for neural networks is known 



to be self averaging in the limit of an infinite sized graph (|Amit et all 



Provost and Vailed . 



1985 



van Hemmen , 



19821 : 



19831) . This occurs because the dynamics is performed over a huge number of 



single neuron excitations. 



1985) in that the 



In practice, our model differs from the neural network model of (|Amit et al.l 
neurons of our model change only from "off" to "on", and cannot flip in both directions as the 
equivalent "spins" do. We therefore tested numerically the self-averaging property of the graph, 
and found that it indeed exists for the random graphs we considered i.e., the fluctuations between 
specific different realizations of the graph are negligible (see Fig. @]below). 



4. The model describes initial growth of activity 



The possibility of turning a neuron "off" is not incorporated in the model because we only 
consider short times. The whole process described by the simulation occurs over a very brief 
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period of time, and therefore a firing neuron keeps its effect on other neurons during the whole 
process. To be concrete, the unit of time in the model and in the numerical simulations is the firing 
of one spike, equivalent in the experiment to about one millisecond. The simulation extends to 
about 50 units, i.e. describes a process that occurs typically for 50 ms. 

In our model a neuron has no internal structure, so that whether it is "on" or "off" impacts 
only the neurons that are its neighbors. The relevant issue is therefore - how long is the effect of 
a neuron's activity felt by its 'typical' neighbor. The experimental facts are that a neuron fires on 
average 4-5 spikes per burst, eac h lasting a millisecond, with about 4-5 milliseconds between the 
spikes ( (IJacobi and Mosesl 120071) ) so that its total active time spans typically 20 milliseconds. 

The post-synaptic neuron retains the input from these spikes over a time scale set by the mem- 
brane decay constant, which is on the order of 20-30 milliseconds. Therefore, after a firing period 
of about 20 milliseconds, there is a retention period of comparable duration. We can conclude that 
the effect of a neuron that has fired can be felt by its neighbors for the total build-up time of the 
burst, about 50 milliseconds. We therefore describe by "on" the long term, averaged effect of the 
neuron once it has begun firing. One caveat to this is that the strength of that effect may vary with 
time, and such an effect is not described within the model. 

We also assume, for simplicity, that all the neurons are available and can participate in the burst 
(no refractive neurons). In the experiment this is equivalent to looking at those bursts that have 
quiescent periods before them, which is often the case. 



5. The role of inhibitory neurons 

In this model all neurons are excitatory; within the "0" or "1" structure of the model, the 
contribution of an inhibitory neuron would be "-1". Thus adding inhibitory neurons amounts 
to increasing m, the number of inputs that must fire before a neuron will fire. This is a small 
accommodation of the model, and does not change the dynamics of burst initiation. 

Below we also model the dynamic s of burst activation observed in the experiments of Eytan 



and Marom (lEytan and Marom , 



20061) . which measured an exponential recruitment at the initial 
stage of the burst. These experiments clearly show that the dynamics are essentially the same for 
cultures both with and without inhibition. Both display an initial exponential growth followed by 
a super-exponent. One small difference lies in the value of the growth, which is larger for dis- 
inhibited than for untreated cultures. But the main difference is seen only after the burst reaches 
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its peak, in the decay of the burst. 



6. The role of noise in burst initiation 



We assume the existence of spontaneous sporadic ac tivity of single neurons in t he cul ture 



In principle, this can be treated as a background noise ([Alvarez -Lacalle and Moses , 



2009). In 



our case we require that a minimal amount / of the culture spontaneously fires, and we look at 
the ability of this fraction of initially firing neurons to initiate a burst. It is possible to initiate 
t he activity with an external voltage V, u sing bath electrodes, as we reported in previous work 



( (|Breskin et al. 



2006; 



Soriano et al 



2008|)). In that case /(V) is determined by the percentage of 



neurons that are sensitive enough to fire at a voltage V. 



HI. RESULTS 



A. First to fire neurons lead bursts and have large input degree 

We use the simulation for an initial look at the recruitment process and to identify those neurons 
that fire first. We use a Gaussian distribution to describe the experimental situation as closely as 
possible, and put the system near criticality, i.e., with / barely above /*, to observe a large range 
of changes in activity. Figure \T\ shows the degree values k of the neurons as a function of the time 
step at which they first fire. It is evident that the neurons that fire first are either the ones that were 
ignited externally or those with high k. This is verified in the lower part of Figured] in which we 
plot, for each neuron, its time of ignition as a function of its in-degree k. It is obvious that the high 
k neurons totally dominate the initial stages of the activity. 

Further information on the distribution of ignition times for different neurons with different 
in-degree k is given in the colored map format of Figure |2l The huge majority of neurons has a 
low k and ignites very late in the burst. The first-to-fire neurons, or Leaders, are few and have 
a wider distribution of in-degree k at a given time step t. The distribution sharpens as the burst 
advances in time. 

To understand this from the model, we concentrate on the neurons within a given A;-class, i.e., 
the group of neurons with k inputs, and examine the probability of a neuron within this group to 
fire, If the average number of inputs k and the threshold value m are large numbers, then we 
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FIG. 1 Top: Average k of all neurons ignited at each time step t from one particular realization of the 
simulation. Bottom: Times of ignition for all 500, 000 individual neurons. It is evident that highly connected 
neurons are the first to fire. For clarity we plot time only from t = 1 and do not show the ~ 1000 neurons 
that were initially ignited at t = 0. We used here / = 0.0033 and p& which for k < fc ta ii is a Gaussian 
centered on k = 75 with width a = 31 and is a power law ~ /c -2 for > k tSL n. We checked that a 
simple Gaussian p^. gives the same qualitative results. 



can neglect the width of the binomial distribution and approximate the error function ^k( m , by 
its limit 0(fc — m/$). This simplifies the dynamics Eq.© into 

/ m \ If for $(t) < m/k 

^(t + i) = / + (i-/)e U-^, = , (5) 

V 1 for$(t)>m/fc 
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time 

FIG. 2 The logarithmic color coding of the number of neurons with in-degree k that ignited at time step t. 
The data are the same as in Figure [TJ 

meaning that under this approximation the whole A;-class fires once $(£) exceeds m/k. In other 
words, the A;-classes are ignited in steps where in each step the classes whose k is in the range 
m/$(t) < k < m/$(t — 1) are ignited. Obviously, the first nodes to be ignited are those with the 
high k. By summation, one finds the iterative equation 

Ht+i) = / + (i-/) £ Pfc e(*-^) =/+(W) |> = /+(W)P (j^j , (6) 

where P(k) = Vk is the cumulative distribution. The dynamics of the approach to the 
fixed point can be graphically described as the iterations between the curves $ and / + (1 — 
/)P(m/$(t)). 

The time continuous version of the iteration equation is 

9(t) = f- ^(t) + (1 - fMm, <&(*)) * / - Hi) + (1 - f)P f^j , (7) 

which can be integrated, at least numerically, to obtain the dynamics of the system. Within this 
approximation, the fc-class firing is given by the step-function ©. For simple collectivity func- 
tions, ^, and simple degree distribution pk © can be integrated analytically. In more complicated 
cases, an iterative scheme, as detailed in Section III.C below, is needed. 

Since these neurons are highly connected, they will statistically be connected to each other as 
well. For the experimentally relevant case of a Gaussian distribution peaked at 78 connections 
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with a width of 25 ( Breskin et ah . 



2006; 



Soriano et al. . 



2008b . more than 10% of the neurons have 



over 100 connections, while about 1% of the neurons have 120 connections or more. Since high- 
k neurons have more inputs and hence a larger probability to receive inputs from other high- A; 
neurons, these highly connected neurons form an interconnected subgraph. We summarize our 
understanding by stating that Leaders are highly interconnected, homogeneously distributed and 
form a sparse sub-network. 

In the Multi Electrode Array experiment about 60 neurons were monitored, and a burst was 
observed to begin with one or two of these neurons. From these initial sites the activity spread. 
Identifying these neurons as Leaders, we reach the conclusion that in every experimentally ac- 
cessible patch of the network that we monitor there is a small number of neurons that lead the 
other neurons in activation. We therefore deduce from the theory that they are part of this highly 
interconnected, sparse sub-network. In the initial pre-burst period nearby neurons are recruited by 
inputs from the Leaders, while in the burst itself all the neurons fire together. During the pre-burst 
a spatial correlation to the Leader exists in its near vicinity, which vanishes as the activity transits 
to the burst. 

We remark here that within our model a node that fires early is highly connected. However, the 
number of connections k and the threshold for firing m are two factors with the same e ffect, and 
they co uld in fact interplay to cause a more complex behavior than we are describing ( Zbindeni 



201 Obi) . One alternative model could hold the number of connections fixed for all neurons, and 
only allow a variation in the number of inputs needed for a neuron to fire m. This would clearly 
bring about a variety of response times of neurons, and could create a subgroup of nodes that 
fire early. If we allow a few neurons to have a low threshold m then those neurons will qualify 
as our Leaders. While there is no evidence to point to a wide variability in the threshold of the 
neurons, there are clear arguments why some neurons may change their threshold in response to the 
activity of other neurons, either reducing the sensitivity (adaptation) or increasing it (facilitation). 
In Section ITTTCl we discuss such a possibility, giving a demonstration of how such a scenario could 
evolve. 
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B. Deducing the connection distribution from the initial growth rates 



1. The experimentally relevant case of an exponential pre-burst 



If the firing order of the neurons is determined by their connectivity, then by observing and 
analyzing the evolution of the burst we may learn about the connectivity of the neurons. We focus 
on the experimental fact that the growth rate of the very first firing is exponential, which leads us in 
the next sections to analyze a particular form of the degree distribution that can lead to exponential 
growth dynamics. 

Our observation that the initial growth of the burst is totally dependent on nodes at the very 
high-k side of the degree distribution gives an opportunit y to find the origin of the initial expo- 



nential growth A(t) 



observed by Eytan and Marom (lEytan and Marorrl 



appears at the very beginning of the burst (i.e., at the pre-burst defined in (|Eckmann et al. 



2006 ). Thi s regim e 



2008)), 



and ends when the maj ority of the network begins to be active and the actual burst (also defined in 



(Eckmann et al. 



20081) ) occurs. During this period the amplitude of activity A{t) grows by a factor 



of about 30, and the value of a is about 0.04 — .05 kHz (a? depends on the time step chosen, which 
is taken to be a millisecond in the experiment (|Eytan and Maroml |2006|) ). After the exponential 
regime comes a phase of faster growth rate, during which the amplitude increases by another facto r 



of about 10. The errors on these factors, taken from Eytan and Marom (|Eytan and Maroml 



2006). 



are estimated to be no more than 10%. We not e that the same exponent ial growth rate is observed 



2007) 



in the experimental data of Jacobi and Moses (|Jacobi and Moses . 

If is known then one can, in principle, extract the in-degree distribution p k , since in the 
random graph scenario nodes with higher k ignite the next lower level, of k— 1-nodes. In particular, 
as we shall now show, an exponential growth rate is obtained for the power law distribution p k = 
Bk~ 2 . We begin by plugging into the approximate dynamics © an exponential time dependence 
$(£) = fe at , where / is the initially lit fraction, and get: 



6 = a$ = /-$ + (l-/)P 



m 



and therefore 



Introducing q = 



m 



Mt) 

P(q) = 



a) • *(t) - / 



/ 



a 



f 



f 



(8) 



(9) 



(10) 
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and we end up with 



Pk 



d 
dq 



P(q) 



m ■ 



q=k 



1 + a 1 



(11) 



This sets the value for B = m- j^j in terms of the growth rate a. We find empirically below that 
the data are best reproduced for a = 0.04, impressively close to the experimental value a = 0.045. 
This indicates that the time steps used in the simulation and in the experiment are similar i.e., the 
firing time of a neuron (simulation time) is very close to a millisecond (experimental time). 



2. The full degree distribution 



However, the distribution obtained above would give an exponential growth at all times until 
the whole network is ignited and $(£) = 1. That is n ot the experimental si t uation. In the data 



2008 



Eytan and Marom. 



of Ey tan and Marom and in that of Jacobi and Moses (|Eckmann et all 
20061) the exponential regime includes a small fraction of the nodes (about 10%). It is followed 
by a faster growth rate, during w hich the remaining nodes fire. We f urther more have measured 



with the percolation experiments (IBreskin et al. 



2006; 



Soriano et al. 



20081) that the distribution 
in a typical culture is well described by a Gaussian, centered on k cen ~ 78 and with a width 
of about a = 25. The average connectivity, as measured by the mean of the Gaussian, was 



20081) . We note that 



shown to increase with the density of plating of the neurons ((Soriano et all 
these percolation experiments measure the fixed point of the firing dynamics and are therefore 
insensitive to any fat tail of the degree distribution, which governs the pre-burst. 

This leads us to the following tailored solution, which combines both these experimental in- 
puts and solves the growth rate problem. We keep the Gaussian distribution for pk over a large 
proportion of the nodes. We have some intuition on why the input degree distribution of nodes 
should be Gaussian: it is essentially determined by the area of the dendritic tree times the density 
of the axons that cross that area (see Discussion). Both the area and the density are expected to be 
random variables in a culture grown on a dish. These random values form a Gaussian distribution 
with a mean and variance that are set by biological processes. 

For simplicity and conformity with the experimental situation, we also demand that no node 
has in-degree less than a minimal k min > m, where m is the number of inputs that need to fire for 
a node to be excited. We therefore begin with the following distribution for small k (A; m i n < k < 
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^tail) : 

(k k cer 

~2a 



p k ~ exp ( J en ) . (12) 



At high k we need to change to a power law distribution pk = Bk~ 2 , and we do this from a 
degree fc ta ii. The value of fctaii, among other parameters, is determined by external considerations 
along with consistency constraints, as detailed below. 



3. Quantitative comparison of model and experiment - setting the parameters 

An impressive experimental fact is the large dynamic range observed in the amplitude of the 
burst, about two and a half decades in total. In the experiment, the amplitude grows in the expo- 
nential pre-burst phase by a factor of about 30, and in the burst itself by a further factor of about 
10. 

In the experiment the large dynamic range and high precision are obtained by averaging over a 
large number of bursts, and can be reproduced in the simulation only if there is a very large range 
of available k values, or else the cascade during which successive A;-nodes ignite each other does 
not last for long enough. To be concrete, we find that we need $ to start from about / = 0.003 
in order to see the amplitude increase by a factor of 300 in total. For the growth to be extended 
in time and to allow sufficient resolution in the simulation, we demand also to be near criticality, 
i.e., f ~ /*. This slows the process by adding only a few neurons that ignite at every time step. 
To obtain such a very long growth time at any other point, away from criticality, would require 
a larger range of k, so by staying near criticality we are actually limiting the range of k to the 
minimum necessary to reconstruct the experimentally determined dynamics. 

To ensure that during the exponential regime $ increases by a factor of 30 while during the 
faster growth it grows by a factor of 10, the transition from exponential growth to the faster, full 
blown firing of the network is designed to occur at $ = 0.1 and / is set at 0.0033. 

Since the entry of the k degree node occurs at a $ = $fc ~ m/k, at the transition from pre- 
burst to burst we have fc ta ii — m/$. Inserting $ = 0.1 and m — 15 gives a characteristic value of 
fc ta ii — 150. This is a considerable distance from the peak of the Gaussian, so that it is justified to 
describe the majority of the nodes by the Gaussian distribution. 

The highest cutoff of the degree distribution is in turn determined by the constraint on the 
integral over the distribution from k ta n to fc max , which should yield a total fraction of 0.1, since that 
is the part of the network that will ignite in the initial, exponential regime, P(A; ta n) = $taii = 0.1. 
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This condition allows us to normalize the cumulative function 



P{k) = $taii • k l • (13) 



k" 1 — k- 

'Hail ""max 



Plugging this P(k) into © yields 



P ( * = ~S 7"^- = $ tail • FT" , ( 14 



where we used fc ta ii = m/$taii- At A; = fc ta ii Eq.(IT4l) yields 

$tail = ~^— f ■ (15) 

Since we have seen that $ tai i = 0.1, this sets consistency demands on a and on /. For $(0) = /, 
we find after some algebra: 

tax = m ■ —J— ■ ( 16 ) 
Since the experimentally relevant values of both / and of a (measured in the appropriate time-step) 
are known and obey / «C a << 1, the approximate relations are $ ta o — f/ot and /c max ~ rn/ f. 
This then sets the pre-factor of the distribution B — m- ^zj — m - We end up with the probability 
distribution function shown in Figure |3j in which a power law tail from fctaii all the way to k max is 
glued onto a Gaussian curve centered on k cen = 75. 

Figure |4] shows our main result, in which an exponential growth rate is reproduced in a simu- 
lation employing the in-degree distribution of Figure [3l This exponential phase is followed by a 
super-exponential phase, in which the majority of the network ignites. The majority has in-degree 
defined by the Gaussian distribution and therefore they fire practically simultaneously. Since we 
tailored $(£), the growth during the exponential phase is indeed by a factor of 30, similar to the 
experimental one. However, the experimental graphs describe the momentary activity while $ is 
the total fraction of active neurons and these do not turn off. The experiment is therefore probably 
better described by the derivative of $, shown in red. It can be clearly seen that in fact $ and $ 
behave very similarly. 

We can also compare the simulations of the network with the numerical solution of the model. 
For this we use the iterative scheme defined by Eq.© to propagate the activity of the network. 
This is shown (in dashed lines) in Figure 0] The excellent congruence of the simulation and the 
mean-field equation gives verification for the use of our mean field model. The success relies 
on the absence of large deviations and insignificance of fluctuations, which is true in our model 
and experiment, due to the benign behavior of the degree distribution and the large number of 
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FIG. 3 The in-degree distribution, p^, plotted in both linear (main plot) and in log-log (inset) coordinates. 
Parameters used for the Gaussian are: k cen = 75, a = 31, k m \ n = 20 while the power law tail p ~ B ■ k~ 2 
goes from fc ta ii = 150 to fc max = 4, 680 and its pre-factor is B = 15.65. This normalizes the distribution 
to integral 1. The log-log scale in the inset highlights the power-law tail, while the linear scale of the main 
plot accentuates the Gaussian that dominates the majority of the population. 

participating neurons. The only deviation from this agreement is at the initial steps, where only a 
small number of leaders are firing. 

In summary, from the quantitative comparison we find that the model has an exponential initial 
transient if the in-degree distribution is mostly Gaussian, with 10% of the neurons in the power 
law tail, and that the highest k can be in the thousands. 

When comparing these results with the experiment, we should remember that only 60 elec- 
trodes are being monitored. The exponential behavior that is observed over a large dynamic range, 
can be resolved since multiple firings at the same electrode are observed with a resolution better 
than 1ms. In the simulation, in contrast, this is modeled by going to high numbers of neurons, 
each of which can only fire once. 
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FIG. 4 Growth rate of burst for tailored distribution from both numerical simulation using 500, 000 neu- 
rons (solid lines) and iteration of Eq.© (dashed lines). Blue lines show overall firing fraction <I> for each 
time step. The red curves show the numerical derivative, <j>. The initial slowing down (the dip in the 
derivative) is due to a clearly evident "bottleneck" in the simulation, during which the firing almost ceases 
to propagate. The parameters used are / = 0.0033, a = 0.03, and k m - m = round (m(l + a)) = 16, 
k max = round (m(l + a)/f) = 4680. 

C. Excitation-dependent threshold 

At the end of Section IIII.AI we noted that the sensitivity of neurons can be changed either by 
varying the number of their inputs, or by varying their threshold. Up to now, we have assumed 
that the firing threshold in the neural network m is a constant that does not change as the burst 
develops, and varied the degree distribution instead. In this section, we examine the impact of 
keeping the connectivity distribution static, while "loading" the recruitment dynamics onto m by 
making it a dynamical variable that depends on the history of neuronal activation. Since varying 
either parameter (m or p k ) will lead to the same results, in principle one could then have any 
distribution p k of input degree, and compensate by varying m. One would then have to verify 
that the necessary variations in m are biologically reasonable and feasible. In the experiment this 
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happens via the competing processes of adaptation and facilitation. Since adaptation would work 
opposite the trend observed in the experiment, we discuss only the possibility of facilitation. 

Facilitation of activity can occur if neurons that are already excited several times are easier 
to excite at the next time. By synaptic facilitation we mean the property of a synapse increasing 
its transmission efficacy as a result of a series of high frequency spikes. We examine here some 
of these effects by introducing, for the sake of simplicity, a threshold which is a function of the 
average firing state, m(<3>). 

Given the time series of the firing fraction $(£) and a presumed degree distribution pi-, one can 
invert equations © or © to obtain 

m{m)=Ht) . K (m±m^i\ (17) 

where K{P) is the inverse to the cumulative function P(k) = J2k™=lPk' (such an inverse function 
exists since P is monotonic). 

It is particularly interesting to ask if the power law degree distribution p k ~ k~ 2 , supports a 
biologically feasible form of m($), following the initial exponential regime where m should be 
constant. In this case the cumulative function is 

p{k) = *L^jsgE , (18) 



max 



where k min and k max are the limits of the distribution. The inverse function is 



K{P) = p.f^+fc-n + jfc-i • (19) 

1 V^min ' "'max/ ~ "'max 



To further advance we need to model the burst itself, which grows exponentially at first, then 
grows even faster, at a super exponential rate and finally saturates when all the network has fired. 
This behavior can be described by the function 



1 — p~ at * 



with t* a parameter to be determined from comparison to the experiment. This kind of burst 
function starts as an exponential $(t) ~ e at and begins to diverge after time steps, but reaches 
= 1 slightly before fully diverging, at t x = U- a' 1 log[l + f(e at * - 1)] ~ t» - (f/a) ■ e at * . 
For this profile $ = a$[l + $/(/(e at * - 1))]. Plugging into CEl]) yields 

m <*w> - TTwrmkn^)] • <21) 
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We see that m($) starts from m(f) ~ A; min /(1 + a) and ends at m(l) ~ k min /[l + a/ {f{e at * — 
1))]. Therefore, m decreases by a ratio of m(l)/m(f) ~ [1 + a/(f(e at * — 1))] _1 , which for the 
experimental parameters is around 5, i.e., from m(f) ~ 15 to m(l) ~ 3, a biologically reasonable 
variation. The actual value of m that we use in the simulation is that of the nearest integer obtained 
by rounding Eq.(fT7T). Figure |5] shows the behavior of the burst as a function of time for the power 
law distribution ~ k~ 2 with variable m($). 




5 10 15 20 25 30 



time step 

FIG. 5 Growth rate of burst for k~ 2 power law distribution and variable m($(t)) using Eq.dTTb. The curve 
is calculated from numerical solution of the iteration Eq.©. Blue line shows overall firing fraction <1> for 
each time step. The red curve shows the numerical derivative, $. The parameter i* used is t* = 40 time 
steps, and all other parameters are as in Figs. [3]and|4] Inset: The threshold m($(t)) decreases during the 
simulation as increases according to d2TT >. The value of m used in the simulation is the integer part of (12)1 
hence the discrete jumps in its value. Until about t = 18 the value m is unchanged and the $(f) and <E>(t) 
profiles are exponential. Then m starts to decrease sharply and induces super-exponential growth of $»(£). 
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D. Summary of Results 



In summary, we have shown here that the experimental situation of an exponential transient 
followed by super-exponential growth can be well described in our model of an in-degree distri- 
bution that is k~ 2 at high k but is Gaussian for the majority of the neurons. The initial transient of 
an exponential is determined directly by the power law tail of the degree distribution of Leaders. 

The A;-degree values needed to describe the data reach a maximum value that is many tens of 
standard deviations from the mean. Although the average value of the degree distribution remains 
in the region of 100, a few neurons (in a network of a million nodes) can have thousands of 
connections. 

We remain with the question of how significant is the need for an exponent —2 in the power 
law distribution, and whether small deviations will change the exponential growth rate. Is there 
any logical or biological reason for this power law to be built up? 

On the experimental side, the search for a few highly connected neurons would be needed. One 
possibility is that Leaders are neurons of a different species then that of the majority. Identifying 
them, investigating their properties and potentially intervening by disrupting their function are all 
important experimental goals. 



IV. DISCUSSION 



A. First-to-fire neurons and Leaders 



In (|Eckmann et all 120081 : IZbindenl l2010allbh Leaders were defined through an intricate math- 
ematical procedure. In particular, this definition allowed for exactly one Leader per burst, which 
ignites a pre-burst, and then a burst. In the present paper, a simpler definition is used, which 
amounts to take into account basically all neurons which fire at the beginning of activity right after 
the initial fraction /. Since in the initial period of the burst there are only very few neurons active, 
the development of the burst depends critically on those neurons. 

Within the QP model, high k neurons activate the low k neurons. So the highest k neurons are 
the ones that trigger the burst. It follows that some of the highest k (in-degree) neurons are both 
first to fire and Leaders. 

Looking only at in-degree is only part of picture. Indeed, high A;-classes are ignited first. How- 
ever, their contribution to the firing propagation depends on the out-degree. Nodes with no outputs 
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may fire early but contribute nothing to the ignition of others. Nevertheless, since we assumed that 
the in- and out-degree are uncorrelated, Leaders are among the early igniting nodes. 

B. How can we get a distribution of in-degree which is Gaussian with a k~ 2 tail? 

An interesting question is what kind of growth and development process would lead to a distri- 
bution of in-degree that is essentially a Gaussian centered on a value of about k ccn = 75, but has a 
tail that goes like k in ~ k~ 2 , and can reach in-degree in the thousands, k max ~ 3, 500. 

We propose the following simplified and intuitive geometrical picture for how k- m and k ont are 
determined. Each neuron in the culture has a spatial extent that is accessed by its dendrites (the 
"dendritic tree") and characterized by a length scale i. The dendrites have no a-priori preferred 
direction, and the dendritic tree is typically isotropic and characterized by a length scale r. Ax- 
ons, on the other hand, go off in one direction, and their length determines the number of output 
connections the neuron will have. The dendritic tree is "presented" to axons of other neurons. If 
the axon of a neuron happens to cross the dendritic tree of another neuron then, with some fixed 
probability (which we take for simplicity to be unity), a connection is made between the two neu- 
rons. The number of in-connections is therefore related to the size of the dendritic tree and to 
the number of axons crossing it, i.e., the density of axons. The number of out-connections of a 
neuron is determined by the length of its axon, the size of the dendritic tree of other neurons and 
the density of neurons. 

There are two corresponding length distributions p(£) and p(r) and a density n that determine 
the number of connections. p(£) and p(r) are the probability distribution of the axon and dendrite 
lengths respectively, while n is the density of neurons per unit area. 

The number of in-connections of a neuron is obtained by calculating the probability of an axon 
emitted from another neuron located a distance t away to cross its dendritic tree. To get the number 
of connections k- m for a neuron with dendritic tree of size r we look for the axons that will cross 
one of its dendrites: 

/•oo 9 r r°° 

k in (r) =n d£27i£—P(£) = 4nnr / d£P{£) . (22) 
Jo " Jo 

Here P{£) = p(£')d£' is the cumulative sum of probability that the length of an axon ex- 
ceeds £ (since it would then cross the dendritic tree). We ignore the slight r dependence of the 
lower limit of the integral, n is the density of neurons per unit area (about 500 neurons per mm 2 ), 
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FIG. 6 Schematic picture of the relation between axon and dendrite lengths i and r to the number of 
connections k m and k out . (a) Two lengths characterize the connections of each neuron: its axonal length £ 
and the typical size of its dendritic tree r. While the dendritic tree is in general expected to be homogenous, 
it can have dendrites that go off much farther than the others, (b) A connection from Neuron A to Neuron 
B will be made if the trajectory of the axon extending from Neuron A will intersect the dendritic tree of 
Neuron B. The probability for that to happen depends on the probability P{£) for A to have an axon of 
length longer than I and on the probability p(r) for B to present a dendritic tree of cross section r. 



and — the angle extended by the dendritic tree as seen from the axon's neuron of origin. 
We now insert for p(£) the Gaussian with power law tail: 

p(£) = A ■ e 2 CT 2 if £ < £ tgil an( i p(f) — £> . £-2 otherwise, with A 3> B normalization 
factors. 

For £ < £ tai i we get 

/•^tail ll'-l f f l ™x 

P{£) = A e — + B dl' = const. - A ■ erf {£) , (23) 

h Je tan 

while for £ tail < I < £ max : 

P(£) = B / £'~ 2 d£' = B[- . (24) 
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The integral over P(£) gives one constant term and one that goes as log(£ max ). Since the maximal 
length is determined by the size L of the culture dish, we remain with a term of log(L). 
We get that the number of in-connections is 

k- m (r) ~ n ■ r • log(L) . (25) 

We are now in a position to ask where the tail of high connections arises. In principle, it could 
arise from fluctuations in the density n. The neural density is theoretically determined by throwing 
down a random set of points on the plane. In the experiment, the neuronal cultures do not exhibit 
clustering to such a degree that would induce so strong a fluctuation with such high density. To 
get the necessary range of a factor of 10 — 30 in k values that the theoretical explanation of the 
experimental data point to, we would need the density to change in a similar manner. This seems 
unrealistic. A further strong argument against fluctuations in density is that these would lead, in 
our picture, to a high number of input connections in one single neighborhood. In particular, this 
would lead to many more high-/;; neurons than the distribution allows for. 

Thus we are led to the conclusion that the power law tail of k in has its origin in the distribution 
of dendritic trees p(r). While the typical dendritic tree is probably a circle of with radius r = 
100 — 200 micrometer, it may have outgrowths in one or more directions that reach as far as 
1 = 1 — 2 mm but can, with small probability, go as far as the dish size L. 

C. Trading off static connectivity distribution for dynamical threshold 

While the possibility of exchanging between elementary dynamics and connectivity statics does 
not come as a surprise, the lesson we take from the results of Section lTlI.Cl is important indeed: One 
cannot distinguish, by observing network dynamics in and by itself, between a static connectivity- 
based mechanism and a mechanism that employs dynamics at the elementary level. While we 
tend to believe it is the connectivity that is dominant, one cannot rule out the neuron's internal 
processes as a possible explanation for the dynamics. 

When only a small fraction $ of the network is active, the chances that a given neuron is 
activated at a high frequency are low. Hence, chances of changes in threshold as a result of 
synaptic or membrane dynamics are low. However, as the active fraction $ of the network grows, 
the chances of a neuron to be bombarded at high frequency become higher; we could then get a 
dependence m($), since changes in threshold (i.e. membrane dynamics) or synaptic efficacy (e.g. 
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facilitation) are expected. 

In many studies of biological networks, this ambiguity is somewhat neglected in favor of a 
more static view, largely due to lack of access to elementary level dynamics. In the case of 
neuronal excitability, single element dynamics is experimentally accessible, and the existence of 
dynamical-thresholds are well documented. Our results in Section UlI.CI indicate that this is a suf- 
ficient explanation to the phenomenon of an early exponential recruitment rate followed by faster 
growth process. 



D. The effect of limiting size: One-dimensional cultures 

Initiation of activity in one-dimensional cultures seems to be ve ry different from the Leaders 



scenario. In one-dimensional cultures, we (Feinerman et al 



2007|) have shown that the activity 



originates locally, at well defined "Burst Initiation Zones" (BIZs) that have a limited spatial extent. 
There are usually a small number of such BIZs, typically one or two per centimeter, that operate 
independently of each other. The BIZs are characterized by a high density of excitatory neurons 
and a low density of inhibitory ones. Firing activity that originates in a BIZ will propagate out as 
a wave-like front with a constant velocity, and invade the rest of the culture until all neurons have 
fired. 

One explanation for the difference in behavior of BIZ and Leaders is in the dimensionality. 
However, the basic argument we presented for the number of input connections relates it to the 
multiplication of the area of the dendritic tree by the density of the axons that cross through this 
area. Since both the radius of the dendritic tree and the width of the line are on the order of 100 
micrometer, there should be no difference in the first factor. As for the density of axons, there is 
no direct information, but also no compelling argument why ID cultures should differ in this from 
2D cultures. 

A different possibility, and the one we b elieve to be correct, is just that there are too few neurons 



in the culture (Tl usty and Eckmann , 



20091) . That would impact on any small culture, both 2D and 
ID. Changing the number of neurons has the largest effect on the realization of the tail of the 
probability distribution pk, since high k values have a low probability and will not be obtained. 
This can completely disrupt the form of the degree distribution. In turn, it also affects the value of 
/*, the fraction of initial firing needed for ignition of the full culture. 

The results of simulating of excitation in varying numbers of neurons are given in Table |U 
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/* 


N 


/c max theory 


&max realized 


0.05 


500 


312 


237 


0.022 


5,000 


709 


615 


0.0076 


50,000 


2,053 


1,795 


0.0051 


100,000 


2,836 


2,512 


0.0038 


500,000 


4,680 


4,680 



TABLE I Study of finite size networks: For small N networks, /* needs to be larger and the experimentally 
accessible fc max does not reach the theoretical prediction. This shows that the ignition process is less efficient 
than for large N. 

We immediately see that indeed /* depends strongly on N. A power law fit indicates that 
/* ~ N~ l l 2 . At about N = 200, 000 the curve flattens out, and reaches the theoretical (N = oo) 
value. The reason for this originates in the constraints imposed between /* and fc max , /* ~ tt ^- 

"-max 

In any realization of finite size N, any k with pk < 1/N is very unlikely to be observed. Since 
Pk ~ k~ 2 both &; max and /* are constrained by the N' 1 / 2 . We can conclude that within the Quorum 
Percolation model smaller cultures require a much larger fraction of initial activity to sustain a 
burst. 
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